Exploring the phase structure of lattice QCD with twisted mass quarks'* 
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The phase structure of zero temperature twisted mass lattice QCD is investigated. We find strong metastabilities 
in the plaquette observable when the untwisted quark mass sweeps across zero. 



1. Introduction 

In the process of approaching in lattice QCD 
the physical point, at which the pion mass as- 
sumes its value as measured in experiment, the 
simulation algorithms suffer from a substantial 
slowing down jll2j which restricts present sim- 
ulations to rather high and unphysical values of 
the quark mass. Moreover the usual Wilson-Dirac 
operator develops unphysical small eigenvalues at 
small values of the quark mass which render the 
simulations more difficult and sometimes even im- 
possible. 

An elegant way out may be the use of Wil- 
son twisted mass fermions 3 with the following 
fermionic action for Nf = 2 mass degenerate fla- 
vors of quarks in the so called twisted basis (x) 

S[ X , X, U] = x(D[U] + too + ^7 5 r 3 )x , (1) 

where D[U] is the standard Wilson-Dirac opera- 
tor, mo is the untwisted quark mass parameter, /j 
is the twisted quark mass parameter and r 3 is the 



third Pauli matrix acting in flavor space. In the 
present paper, unless otherwise stated, the lattice 
spacing is set to unity: a = 1. The twisted mass 
H serves as a natural infrared regulator for the 
low lying eigenvalues of the Wilson twisted mass 
operator since 



det (D[U] +to +fiij 5 T 3 ) 
= det((D[U] + m Q )(D[U} 



too) 



V) 



(2) 



*Presented by C. Urbach and F. Farchioni. 



Note that on the l.h.s. the operator is the two 
flavor operator, while on the r.h.s. D[U] + too 
is only the one flavor part. As for the gauge ac- 
tion, the usual Wilson plaquette action is used. 
Note that the bare quark mass too is often rep- 
resented by the hopping parameter k defined as 
k = (2m + 8)- 1 . 

In addition to the infrared cut-off the Wilson 
twisted mass formulation allows to get full O(a) 
improvement for correlation functions and de- 
rived quantities with no need of additional coun- 
terterms [415] . provided too is set to its critical 
value TO cr it and the value of [i is kept constant as 
a -> 0. 
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Figure 1. Thermal cycles in k on 8 3 x 16 lattices 
at /3 = 5.2. The plaquette expectation value is 
shown for: fj, = 0.1 (A); fi = 0.01 (B); ^ = 
(C). The triangles refer to increasing k- values, the 
diamonds to decreasing ones. 



Studying the phase structure of lattice QCD 
should be a pre-requisit before starting to ex- 
tract physical results. Indeed, in simulations of 
Nf — 3 clover improved Wilson fermions || and 
of Nf = 2 non-perturbatively improved Wilson 
fermions [2] metastabilities and hysteresis effects 
were found. In this contribution, or in more detail 
in we discuss results for the phase structure 
of Nf = 2 Wilson twisted mass fermions. 

For our simulations of full QCD we have imple- 
mented two independent algorithms: The Two- 
Step Multi-Boson algorithm (TSMB) and the Hy- 
brid Monte Carlo algorithm (HMC), both with 
even-odd preconditioning. For the HMC we use 
in addition the Hasenbusch trick [SJ. We checked 
that we get the same results with both algo- 
rithms. 

We have performed simulations primarily at 
(3 — 5.2, but we also have some data at /3 = 5.3 
and (3 = 5.4. The lattice sizes are 8 3 x 16, 12 3 x 24 
and 16 3 x 32. 

2. Thermal cycles and metastable states 

We started our investigation of the phase dia- 
gram by performing thermal cycles in k. Fixing 
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Figure 2. Monte Carlo history of metastable 
states at /3 = 5.2, fi = 0.01 and k = 0.1717. In 
the upper plot the lattice size is 16 3 x 32 and in 
the lower plot it is 12 3 x 24. 



the values of — 5.2 and fj, we incremented k 
from a starting value K sta rt until the final value 
Kfinai was reached and then reversed the proce- 
dure. At each value of k 150 configurations were 
produced - without performing further interme- 
diate thermalization sweeps - and averaged over. 

In fig. H we show three such thermal cycles, 
performed at fi = 0, /i = 0.01 and /i = 0.1 from 
bottom to top. In the cycles signs of hysteresis 
effects can be seen for /i = and /.i = 0.01 while 
for the largest value of \i = 0.1 such effects are 
hardly visible. Hysteresis effects in thermal cycles 
may be signs of the existence of a first order phase 
transition. However, they should only be taken as 
first indications. 

Guided by the results from the thermal cycles, 
we next performed simulations at fixed values of 
fi and k, starting with ordered and disordered 
configurations, staying again at (3 = 5.2. 

In fig. |21 we show the Monte Carlo time evo- 
lution of the plaquette expectation value for two 
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Figure 3. Pseudoscalar mass as a function of 
1/(2k) at /3 = 5.2 and fi = 0.01. 



_l I I l_ 



12 3 x 24 low I □ I 
16 3 X 32 low I A | . 
12 3 X 24 high O 
16 3 X 32 high h^H 
_l I I I I 



2.895 2.9 2.995 2.91 2.915 2.92 2.925 2.93 2.935 2.94 2.945 
VP*) 

Figure 4. m^ CAC as a function of l/(2/t) at P 
5.2 and fi = 0.01. 



sets of parameters at P — 5.2. In the upper plot 
the lattice size is 16 3 x 32 and in the lower plot 
the lattice size 12 3 x 24, both with fi = 0.01 and 
k = 0.1717. We find coexisting branches with 
different values of the average plaquette with a 
rather large gap in between. As can be seen the 
gap is not decreasing with increasing lattice size 
and therefore this behavior cannot be ascribed 
to any finite volume effect. Furthermore we ob- 
served this phenomenon in simulations with fi = 
and with both algorithms. We therefore conclude 
that the existence of metastable states is a generic 
feature of lattice QCD in this formulation. 

3. Pseudoscalar and quark masses 

In order to study the physical properties in the 
two metastable states we measured the (charged) 
pseudoscalar meson mass and the untwisted 
PCAC quark mass. We obtained the pseudoscalar 
mass from the pseudoscalar correlation function 
in the %-basis while we measured the untwisted 
PCAC quark mass from the axialvector current 
in the ^-basis: 

m p C Ac = (&Hxw5 I rx(x) Q T (y)) (3] 

™ X 2(x^ 7 5X(z) 0*(y)) 

Here T is a suitable operator that we have 
chosen to be the pseudoscalar density O t = 
X 2 2~75x(x), d* is the lattice backward deriva- 
tive defined as usual and r ± = t\ ± iT2- One 
can show that in the continuum limit a — » the 



quantity is asymptotically proportional 

tO TO - TOcrit- 

In fig. |3 we show the pseudoscalar mass in lat- 
tice units as function of 1/(2k). We observe that 
the "pion" mass is rather large and the most strik- 
ing effect in the graph is that it can have two 
different values at the same k. If we consider 
the quark mass to^ cac in fig. 0] we see that in 
the phase with a low plaquette expectation value 
the mass is positive while for high values of the 
plaquette expectation it is negative. Values of 
m PCAC with opposite sign coexist for some val- 
ues of k. From this fact one can safely argue that 
the value of 1/ (2/c cr i t ) lies between 2.91 and 2.92. 

Figures 121 II clearly reveal that for fi = 0.01 
metastabilities show up in the quantities we have 
investigated, such as mps, m^ CAC and the aver- 
age plaquette, if mo is close to its critical value. 

4. Determination of the twist angle 

The knowledge of the twist angle to, the polar 
angle in the plane of the untwisted and twisted 
mass, is important e.g. when comparing lattice 
data with analytical predictions [H] from Wilson 
chiral perturbation theory (WChPT). We present 
here a method which allows to determine the 
twist angle only on the basis of symmetry. 

Following P|, we introduce 2 the twist angle uj 
as the chiral rotation angle between the renor- 
malized (physical) chiral currents , A° and the 
bare bilinears of the x-helds V®, with renor- 

2 In this definition of the twist angle was called a. 
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malization constants Zy and Z a- Thus we have 

V" = Z V V« cos u + e ab Z A A\ sin u (4) 
Al = Z A Al cos uj + e ab Z V V* sin a; , (5) 

where only charged currents are considered 
(a=l,2). The twist angle u> is related to the ratio 
of the renormalized twisted and untwisted masses 
entering the chiral Ward identities 0: 

w = arctan (/ifl/mfl) . (6) 

We define in addition the two auxiliary angles 

uiy — arctan(Z J 4^y 1 tan uj) , (7) 
uja = arctan(ZyZ_^ 1 tanw) . (8) 

In terms of ujy, uja eqs. (@J), © are written 

V£ = My (cos + e ab sin uj v Al) (9) 

Al = M A (cobwa^J + e afc sinw A V; b ) , (10) 

where the overall multiplicative renormalization 
reads 



A'' 



Z x 



x 



(11) 



cos u>x VI + tan uiy tan uja 

From 0, (BJ it follows: 

a; = arctan (Vtanwy tanw^) . (12) 

uiy and wa can be directly determined by im- 
posing parity conservation for suitable matrix el- 
ements; e.g., with P ± (x) — X I 2~75X( a; ) : 

^{A+(x)V -(y)} = Y / {V+(x)P-(y))=0. (13) 



These equations admit the solution 

tan a; a = 

^ S ,g{A+{x)V -{y))+t^ v ^{At{x)A^{y)) 
E*y(V + (x)V - (y))^nu> v j: s JV + (x)A^ (y)) 



tanwy 



-iE s ,y(v + (x)P-(y)) 



(14) 



Alternatively, a determination of ujy and uja 
is given by the vector and axialvector Ward 
identities, respectively; e.g. in the vector 



case, by enforcing the Ward identity with the 
insertion of some appropriate operator 0{x): 

(d;v+(x)6-( y ))=o. 

Once ujy and uja are determined, the twist an- 
gle uj is obtained by eq. l(T2"|) . The method de- 
scribed above for determining the twist angle can 
also be used in case of simulations with partially 
quenched twisted mass quarks. The estimate of 
uj is of course affected by 0(a) ambiguities. For 
^cx/irt^O, |w |=7r/2 signals m =m cri t. 

We determined the twist angle for sets of con- 
figurations at /3—5.2, fi— 0.01 in the positive and 
negative quark mass branches; results are re- 
ported in fig. El 

Owing to eq. 10 uj should approach the value 
uj = tt/2 from above or below, respectively, when 
the untwisted quark mass m^ CAC oc to_r gets close 
to zero from negative or positive values and fi oc 
/ir is kept fixed to a nonzero value. For [i = 0.01, 
outside the region of metastabilities, we observe a 
trend of u> consistent with the above expectation 
(see fig.|SJ|. In the metastability region (2n)^ 1 = 
2.91 — 2.92 the measured values of u lie above 
or below tt/2 in the high or low plaquette phase, 
respectively, reflecting the behavior of m^ CAC in 
fig. The values of uj are far away from tt/2. 

Eqs. l(§|l- i|ll|) allow to determine the physical 
currents and A° up to the usual multiplica- 
tive renormalizations Zy and Za- The physical 
PCAC quark mass was computed from the Ward 
identity for the physical axialvector current, anal- 
ogously to © 



PCAC _ N^M A 



{ d;z^A±(x) 6*(y)) 

2(P±(x) OT(y)) 



(15) 



We determined the pion decay constant from 
the physical axialvector current. The results are 
plotted in fig. versus the physical PCAC quark 



5. Interpretation of the results 

The observed metastabilities characterizing co- 
existing phases with opposite sign of the un- 
twisted quark mass support the picture of a first 
order phase transition (PT) in the niQ-fj, plane for 
small values of [i and mo close to m cr i t . Also the 
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Figure 5. Twist angle as a function of l/(2re) at 
/3 = 5.2 and /i = 0.01 in the two branches. 
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Figure 6. Pseudoscalar decay constant as a func- 
tion of the physical PCAC quark mass at j3 = 5.2 
and \i = 0.01 in the two branches. 



observed gap in the average plaquette finds a nat- 
ural explanation in the presence of a first order 
PT, when the effects of the spontaneous breaking 
of chirality are considered in combination with 
the explicit breaking of the symmetry due to the 
Wilson term in the fermionic action. 

A first order PT for Wilson fermions at small 
0(a 2 ) quark masses is one of the two possi- 
ble scenarios predicted by chiral perturbation 
theory with inclusion of lattice corrections [T2*] 
(WChPT). The other one is the well-known Aoki 
phase, known to be realized at strong gauge cou- 
plings (see 531)- Which of the two scenarios ap- 
plies depends on the sign of the second order 
(0(a 2 )) corrections in the effective potential. The 
generalization of the results of ref. to the case 
of twisted mass fermions was recently worked out 
independently in and the first order 

PT extends to \i ^ getting weaker and weaker 
for increasing li; at the endpoint li c = 0{a 2 ) the 
PT becomes a critical point. 

The metastability of two phases with positive 
and negative untwisted quark mass has an in- 
terpretation, at a microscopic level, in terms of 
the properties of the eigenvalue spectrum of the 
Wilson-fermion matrix; studies in this direction 
are in progress JS] (see ^H] f° r a numerical study 



of the case /x=0). The qualitative picture ^B] is 
that the tunneling between the two phases im- 
plies a massive rearrangement of small eigenval- 
ues. This process is suppressed by the zero of the 
fcrmion determinant when the eigenvalues move 
close to the origin in the complex plane. For lij^O 
and large enough, the metastability is lifted since 
the condition |ImA| > /i implies the depletion 
of the spectrum around the origin. This conclu- 
sion is in agreement with the prediction of twisted 
mass WChPT (see above). 

6. Conclusion and outlook 

We investigated the zero temperature phase di- 
agram of twisted mass Wilson fermions. The 
inclusion of a twisted mass term made it easier 
to explore the region of small (untwisted) quark 
masses. Metastabilities, signaled by different val- 
ues of the plaquette, were visible in this region in 
thermal cycles and long-living metastable states. 
The two metastable branches are characterized 
by opposite values of the untwisted quark mass. 
The observed phase transition is rather strong: 
instabilities start already for heavy (~740 MeV) 
pion masses. This is also reflected in the values 
of the twist angle in the two metastable branches, 
which are far from cj = tt/2. 
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We interpret our findings in terms of the ex- 
istence of a first order phase transition line in 
lattice QCD with Wilson fermions; this scenario 
is supported by Wilson chiral perturbation the- 
ory. We argue that a phase transition for small 
quark masses is accompanied at a microscopic 
level with the properties of the spectrum of the 
Wilson-fermion matrix in proximity of the origin 
in the complex plane. 

The existence of a first order phase transiton 
line for small quark masses poses the question, 
how a critical bare mass can be estimated in prac- 
tice. One possible choice is to identify m C rit with 
the value of mo where the pion mass is minimal. 
With the untwisted mass set to this value and for 
// > He — 0(a 2 ) the vacuum is unique. 

There are still many aspects which must be 
clarified in order to get a more detailed picture. 
First, how the present metastabilites are related 
with those observed in the literature with differ- 
ent kind of actions and number of flavors |17I6I2| 
(see also ^H] for a recent discussion). We shall 
start by investigating the influence of the lattice 
gauge action on the phase structure near zero 
quark mass |15|. Preliminary results at lattice 
spacing a ~ 0.2 fm show that replacing the Wil- 
son plaquette action by the DBW2 action makes 
the minimal pion mass and the jump of the aver- 
age plaquette substantially smaller. 

Another direction of investigation is the depen- 
dence of the metastability on the gauge coupling 
/3. In the continuum limit (/? — ► oo, a — > 0) the 
minimal pion mass and the jump of the average 
plaquette are expected to vanish and the first or- 
der phase transition line is expected to shrink to 
a first order phase transition point. 

Finally, the most important question is 
how phenomenology can be done with Wilson 
fermions in presence of the metastability. In other 
terms, one has to determine the lightest quark 
mass that can be simulated in a stable phase, 
for values of lattice spacings accessible to present 
simulations. 
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